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We present a complete prescription for the numerical calculation of surface Green's functions and 
self-energies of semi-infinite quasi-onedimensional systems. Our work extends the results of Sanvito 
et al. [l|] generating a robust algorithm to be used in conjunction with ab initio electronic structure 
methods. We perform a detailed error analysis of the scheme and find that the highest accuracy 
is found if no inversion of the usually ill conditioned hopping matrix is involved. Even in this case 
however a transformation of the hopping matrix that decreases its condition number is needed in 
order to limit the size of the imaginary part of the wave-vectors. This is done in two different 
ways, either by applying a singular value decomposition and setting a lowest bound for the smallest 
singular value, or by adding a random matrix of small amplitude. By using the first scheme the size 
of the Hamiltonian matrix is reduced, making the computation considerably faster for large systems. 
For most energies the method gives high accuracy, however in the presence of surface states the error 
diverges due to the singularity in the self-energy. A surface state is found at a particular energy if 
the set of solution eigenvectors of the infinite system is linearly dependent. This is then used as a 
criterion to detect surface states, and the error is limited by adding a small imaginary part to the 
energy. 

PACS numbers: 72.10.Bg,73.63.-b,71.15.-m 
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I. INTRODUCTION 



The electronic transport properties of quasi- 
onedimensional (ID) systems, described by a localized 
orbitals basis set, can be calculated using the nonequi- 
librium Green's function (NEGF) method .^i^iii^ The 
system is usually divided into two semi-infinite left- and 
right-hand side leads, and a scattering region joining 
them. The effect of the leads onto the scattering region 
is taken into account by the so called self-energies 
(SE), which can be calculated from the surface Green's 
function (SGF) of the semi-infinite leads. These can be 
obtained either with recursive methods^'ii^'^ or by using 
a semi-analytic formula i^'^'^^i^^'^^ Recursive methods are 
affected by poor convergence for some critical systems, 
typically when the Hamiltonian for the leads is rather 
sparse. Semi-analytical methods instead bypass those 
problems by construction, however major difficulties 
arise if the hopping matrices are singular or, more 
generally, ill conditioned. Unfortunately the condition of 
the Hamiltonian is set by the electronic structure of the 
leads and by the unit cell used, and thus it is largely not 
controllable. For this reason an algorithm that performs 
under the most generic conditions is highly desirable. 
Here we present an improved semi-analytical method 
that overcomes these limitations and thus represents a 
robust algorithm for quantum transport based on ah 
initio electronic structure. 

In the first part of the paper the extended algorithm for 
the calculation of the SE is presented. First the construc- 
tion of the Green's function of an infinite ID system as 
derived in reference [l| is recast into a more general form 
based on the notion of a complex group velocity. Then we 
present an extension of such method to the calculation 
of the SGF and SE. The new algorithm is defined also 



for the case of singular hopping matrices. This largely 
improves the numerical accuracy. However we find that 
even such an improved scheme sometimes fails if the hop- 
ping matrices are close to being singular. We overcome 
this problem by performing a transformation of the hop- 
ping matrix that reduces its condition number k, defined 
as the ratio between its largest to its smallest singular 
valueJ^i^ This transformation limits the maximum ab- 
solute value of the imaginary part of the Bloch wave- 
vectors, increasing both accuracy and stability. Two ap- 
proaches are presented, the first is based on a singular 
value decomposition (SVD), which is also used to signif- 
icantly reduce the dimension of the Hamiltonian, while 
the second consists in adding a random noise matrix. 
This extended scheme is implemented in the NEGF ah 
initio transport code Smeagol^^ based on the density 
functional theory (DFT) code SIESTAJ^ 

In the second part of this work we present three exam- 
ples of calculations performed with our new implemen- 
tation. We compare the results to the ones obtained by 
using the original method of reference finding a con- 
siderable improvement. However, although the algorithm 
appears very robust, our detailed error analysis reveals 
that for a given system the accuracy is lost at some spe- 
cific energies. This is caused by the divergence of one of 
the SE eigenvalues. The physical origin of this behavior 
lies in the presence of surface states very weakly coupled 
to the semi- infinite leads. Surface states appear when- 
ever at a given energy the set of Bloch functions (with 
both real and imaginary wave-vectors) for the infinite 
quasi-lD system is linearly dependent. In the simplest 
case this corresponds to two Bloch functions being equal 
inside the unit cell. A small imaginary part is thus added 
to the energy in a small energy range around the surface 
state. It is shown that this has little effect on the trans- 
port properties in the high transmission regime, whereas 



2 



Hr 



Hr 



Hn 



Hr 



Hn 



FIG. 1: Schematic representation of the system with onsite 
Hamiltonian Ho and hopping Hi . The overlap matrix has the 
same structure. 



for low transmission it has a substantial influence on the 
results. Crucially only a very small imaginary part is 
used, and moreover this is added only around the energy 
of the surface state and thus the error can be carefully 
controlled. 



II. RETARDED GREEN'S FUNCTION FOR AN 
INFINITE SYSTEM 

Following the scheme introduced in reference [l[ the 
construction of the retarded Green's function for an infi- 
nite quasi-lD system is now recalled. This is the starting 
point for the calculation of the SGF. It is assumed that 
the Hamiltonian is written over a localized orbitals basis 
set and that the interaction has finite range. The size 
of the unit cell can be chosen to guarantee interaction 
only to the first nearest neighboring unit cells. The total 
Hamiltonian of the system H^z' (the integers z and z' 
label the unit cells) can then be written as 



Ho 6, 



■ z' + l, 



(1) 



where Hq, Hi and H^i arc NxN matrices, with being 
the number of orbitals comprised in the unit cell (see 
figure[T]). If time-reversal symmetry holds then Hq = Hq, 
and H-i = h\, however the solutions presented here are 
valid also in the more general case when Hq ^ Hq and/or 
7^ hI. We further assume that the overlap matrix 
Szz' has the same structure and range of the Hamiltonian 

Szz' = •S'o Szz' + Si Sz^z'-l + S-i Sz,z' + 1, (2) 

where 5*0 , and S- 1 are again NxN matrices with the 
same meaning of their Hamiltonian counterparts. 



A. Bloch states expansion 

The solutions of the Hamiltonian equation for the 
associated infinite periodic system J^z' ^^z' i^z' = 
E Szz' i^z' are Bloch functions ipz = e^'^^cf), where 
ipz and 4> are A^-dimensional vectors and k is the wave- 
vector, which in general is a complex number. For a 
given real or complex energy E there are 2N solutions 
with wave-vectors fc„ and corresponding wavefunctions 
d>„ . Each of them satisfies 



(Ha+Hie'''" + ff_r 



E{Sn + Sie'''"+ S^i e' 



(3) 



If we define Kq. = H^ — ESa, {a 
above can be rewritten as 



-1,0,1), the equation 



(4) 



where the additional index R denotes explicitly that the 
solution is a right eigenvector. The corresponding left 
eigenvector (/)l,ti satisfies 



4. 



A'_ie 



= 0. 



(5) 



Time-reversal symmetry gives (/)l,„ = 0L(fcn) = 4>R{kn), 
so that in the case of real fc„ (propagating states) left 
and right eigenvectors are equal. For complex k„ left 
and right eigenvectors are diff'ercnt, describing left- and 
right-decaying states. The sets {fc„}, {0r,„} and {(f>L,n} 
that satisfy eqs. ([4]) and ([5]) at a given energy can be 
found by solving a quadratic eigenvalue proble m^^i^^ of 
the form 
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where 



R,n 



L.n 



2 K-l 



(6) 
(7) 

(8) 
(9) 



Here 1 and (D are respectively the NxN unit and zero 
matrices. The normalization constant is the square root 
of the complex group velocity = dE/dkn [h = 1) 
equal to 



9R,n, 
^R,n- 



(10) 
(11) 



In the following we assume that the eigenvectors 0r.„ 
and 4'-L,n are always normalized to give Z„ = 1. If time- 
reversal symmetry holds then w(A:*) = w^, so that if the 
imaginary part of kn is zero the group velocity is real. 
Note that, at variance with reference [l|, eqs. ^ and 
([7|) avoid the inversion of Ki, so that they eliminate a 
possible source of singularities in the calculation of fc„, 

0R,n and (^L,n- 

The full sets of left {$L,n} and right eigenvectors 
{<I>R^„} form a complete and orthogonal basis. The or- 
thogonality relation is 



L,7l 



Ki (D 



R,m 



(D 1 

where c„ is a constant. This leads to 

^R.m — ^r) 



(12) 



(13) 
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For n = m this equation is only satisfied if c„ = 1, in 
which case it corresponds to the definition of With 
the chosen normaHzation the basis is therefore orthonor- 
mal. The corresponding completeness relation then reads 



2N 



t ( Ki 
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and provides the three following useful relations 

27V 



E 



n=l 



n=l 



V, 



n 



I 1 + 



71^1 



(14) 

(15) 
(16) 
(17) 



Note that in eqs. (|14m7[) the sums run over all 2N solu- 
tions. If Ki = kI^-^ and Kq = eqs. and ([T7]) are 
equivalent. 



The inversions in eqs. (pO|) are usually well defined, unless 
Q and Q do not have full rank. We will return on this 
aspect in section IVll for the moment we assume that the 
duals can always be constructed. 

The Green's function calculated in reference [ij is then 

_„/ E^Li'^R,"^''"^^-^'^^,^^-' ^>^' (21) 



. Eli^R,ne'^"(^-^')0^„„y-i z<z', 
with the matrix V = = .g^Q^ given by 




— ikn 




We now introduce the right transfer matrices^i^iii^ Tr and 



N 



Tr 



Tr 



ri=l 
N 



E "^R^' 



"R,n- 



(23) 
(24) 



B. Green's function 

The retarded Green's function g^z' of the system is 
defined by means of the Green's equation 



'qzz, [{E + i5) 



(18) 



with (5 — > 0"*" real. In what follows we present and expand, 
by using left and right Bloch functions, the solution to 
eq. (dH]) given in reference [l| only in terms of the right 
eigenvectors 0r. First we divide the 27V vectors 
into A'' right-going states with either Im(fc„) > (right 
decaying) or Im(fc„) = and w„ > (right propagat- 
ing), and N left-going states with either Im(fc„) < (left 
decaying) or Im(fc„) = and u„ < (left propagating). 
As a matter of notation in order to distinguish left- from 
right-going states, in what follows wc indicate the right- 
going states with k, 4> and v, and the left-going states 
with a bar over these quantities, i.e. k, (j) and v. 

As in reference [l| we introduce the duals cpR^n of the 
right-going states 0r^„ defined by ^|^ „(^R,m = (5„,„, and 
the duals 0r,„ of the left-going states ^r^„ defined by 
^R n'^R,™ ~ ^nm- If WC dcfinc the matrices Q and Q as 



Q = ( 0R,1 (f>R,2 ■ • ■ (t>K,N ) , 
Q = ( 4>R.l 4>R.2 ■ ■ ■ (t>R,N ) , 

then the duals can be obtained by simple inversion 

( ^R,l 0R,2 •■• (f>R,N ) = {Q^^y , 
{ 4>R,1 4>R,2 ■ ■ ■ (t>R,N ) = {Q 



(19) 



(20) 



Note that both Tr and Tr have eigenvalues with complex 
modulus < 1. For an integer z the following relations 
hold 



N 



(Tr/ =5]</.R^„e^'="^4^„, 

n=l 



Tr 



N 
n=l 



(25) 



which allow us to write the Green's function of equation 
^ as 



9z 



(Tr)^ ^ 500 



[ (Tr) 500 
In the same way V is rewritten as 



V = 5o"o' = K-i {T, 



R 



z>z' 
z < z' 



Tr 



(26) 



(27) 



Note that although the matrices Tr and Tr are in gen- 
eral well defined, the inverse of these matrices is not. In 
fact if A'l and K-i are singular there are some fc„ with 
Im(fc„) — > oo, so that e**^" = (see section llV Ap . In this 
case Tr does not have full rank and is therefore singu- 
lar. The same argument holds for Tr. Equation ((27)) can 
therefore be used only if the matrices Ki and K^i are 
not singular. 

A possible way for overcoming such limitation is by 
using an equivalent form for the Green's function based 
on the left and right eigenvectors. The starting point 
is the relation that will allow us to find the con- 
nection between the duals and the left eigenvectors. 
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Eq. ((T5|) contains a sum over both left- and right- 
going states. By moving the contribution of the left- 
going states to the right side of the equation wc obtain 



N 



N 



, — -. — / , — '" = B, where we have 

introduced the auxiliary matrix B. By multiplying B 



from the left with either 0p or 0^ we obtain respectively 
= J-(h] B-^ and $i =^ B-\ The ma- 



trix B is determined by inserting these relations into eq. 

and by using eq. P?|) . the result is i? = goo- The 
relation between the dual basis and the left eigenvectors 
is therefore 



1 t -1 

i (7 J), 



1 



^L5o"o'- (28) 



This result allows us to rewrite the Green's function of 
eq. ([2T|) in a shorter form 



AT 



AT 



z> z' 
z < z' . 



(29) 

This result represents a generalization to complex en- 
ergies and to systems breaking time-reversal symmetry 
of the solution given in references 2Clll2l| for Hermitian 
Hamiltonians, real energy and an orthogonal tight bind- 
ing model. This derivation shows that the Green's func- 
tion can be cquivalcntly expressed by using the right 
eigenvectors and their duals (eq. ([2T|) ). or both the right 
and left eigenvectors (eq. (US])). It is thus possible to 
move from one representation to the other through eq. 
([28|) that relates the duals to the left eigenvectors. One 
can then decide which representation to use, depending 
on the specific problem investigated. We note that eq. 
((29)) has the benefit that goo can be calculated also in the 
case where the two matrices and K_\ arc singular. 
For those fc„ where Im(fc„) oo the group velocity bc- 



The structure of eq. ((33|) is the same as that of eq. (|26p . 
with the difference that now goo is multiplied to the left 
of the transfer matrix. Finally we extend eq. (^7)) and 
present four equivalent relations for the inverse of goo 

9o,' = K-i {Tn' - Th) = Ki (Tr ' - Tr) 

= {T^'-n)K^, = {f^'-n)K,. (34) 

The second of these relations can be shown by multiply- 
ing eq. (jl6p by from the right and then by using eq. 
(|28p . In the same way the third and fourth equations 
can be obtained by multiplying equations ^TE\\ and (fTT)) 
by from the left. 

In the following we will use mostly the quantities ex- 
pressed in terms of the right eigenvectors only, however 
the same conclusions can be derived using the left eigen- 
vectors. 



C. Density of states 

As an example of the use of the Green's function in 
the form of eq. (j29p we determine the spectral function 
A and the density of states (DOS) of the infinite quasi- 
ID system for the special case where the Hamiltonian 
and the overlap matrices are Hermitian. The spectral 
function is defined as^ 



A-zz' — 



H9-9 



9z 



The DOS pz projected on the unit cell z then is 



1 

:^Tr 
27r 



^ ^ A-zz' Sz'z 



(35) 



(36) 



comes Vn 



'i nKo<pK.n and is therefore well defined ^y using eq. @ this becomes 



{vn = -i 4>i^n^o4>R,n for Im(fc„) -oo). 

As a matter of completeness we show that a repre- 
sentation entirely based on the left Bloch functions and 
their duals (/)L.n and 0L,n is also possible. By multiplying 

eq. respectively by (/)l,ti and (^L,n from the right we 
obtain the two relations 



Pz - — Tr [AzzSo + Az,z-iSi + Az.z+iS^i] . (37) 

ZTT 

In general the main contribution originates from the first 
term in the sum, which can be interpreted as the onsite 
DOS Pz 



^L,7i — -. 9oo'f'R,n 



1 



-IV; 



-5ooVr,"- (30) 



2tt 



Tr [AzzSo] . 



(38) 



Again the left transfer matrices Tl and Tl are defined as 

N 

Tl = E^L,„e^'^-"<_„ (31) 



Tl=l 

N 



(32) 



and the Green's function of eq. ([29| can be rewritten as 



9z 



.900 {TlT ^ 
900 {TlY 



z>z' 
z < z' 



(33) 



We now calculate A and p for the special case where 
K_i — k\ and A'q = Kq. In this case for Im(fc„) = we 
have (/)L,Ti = </'R,Ti, whereas if Im(fc„) ^ then 0L.ri = 
'phikn) — 4>Rikn)- In the same way for Im(fc„) = 
we have (jji^n = 4'R,n, whereas if Im(fc„) 7^ then 
^L,n = 4>L(kn) = (^R(fc*). Therefore for each right de- 
caying state with Im(fc„) > there is a left decaying 
state with fc„ = fc* and v{kn)* = v{kn)- Using these 
relations when inserting the Green's hmction of eq. ([29]) 
in the definition of Azz' , the contribution from all the 
decaying states cancels out. The only remaining contri- 
butions come from the propagating states, also denoted 
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as open channels. For these fc* = kn, kn = — fc„ and 
v{kn) = —v{kn)- With these constraints, and by using 
eq. (|29|) . the spectral function becomes 



Wopcn 

E 



ikn (z — z') 



^ — ikji {z — z ) 



(39) 

where A'opcn is the number of open channels (number of 
Bloch functions at a given energy with real positive k 
vector). If there are no open channels Azz' = and the 
Green's function is Hermitian. Finally, by using eqs. (|39p 
and ([57]) . and the fact that the eigenvectors are normal- 
ized so to give In = 1 (sec cq. the DOS at the site 
z = is simply 



Po = - E — (40) 



This is the well known result for the DOS of infinite pe- 
riodic ID systems.— 



III. SURFACE GREEN'S FUNCTION AND 
SELF-ENERGY 

The retarded Green's function go, for a quasi-periodic 
system, where the left and right sides are separated at 
the position z ~ (the left-hand side part extends from 
z = — ootoz = — 1 and the right-hand side part from 
z ~ 1 to z ~ oo, with no coupling between the cells 
at z = —1 and z = 1), can be constructed from the 
Green's function g for the infinite chain as demonstrated 
in reference [l| 



9zz 



9z0 3oo 9oz'- 



(41) 



The left-hand side SGF is then defined as Gl = gs,~i,~i, 
and the right SGF as Gr, = 5s, n- The SGF can be 
obtained by using eq. 



(42) 



Gl = (l-7RTR)goo, 
Gr = {l-TB.Tn)goo- 



This corresponds to the form derived in reference [l| 
This result can be simplified by using the relations ([34] 
for goo to 



Gl = K^^, 
Gr = Tr KZl- 



(43) 



These equations unfortunately are only defined if Ki and 
-ftT-i are not singular. The same problem however does 
not afi'ect the left and right SE, Sl = K-i Gl Ki and 
Er = Ki Gr K-ir- since they simply are 



Sl 
Er 



(44) 
(45) 



In complete analogy the same expressions obtained by 
using the left transfer matrices are Sl = Ti^Ki and 
Er = TLi^-i. This result is equivalent to those obtained 
in references [sHol flolflllfT^ and derived with different ap- 
proaches, demonstrating the equivalence of those to our 
semi-analytical formula. Since NEGF-based transport 
codes simply require El and Er , our scheme allows the 
calculations of system with arbitrarily complicated elec- 
tronic structure. A schematic tree diagram describing 
the steps involved in obtaining the SE is shown in figure 
[5] (basic algorithm). 

Eqs. and (|^ demonstrate that the SE can be 

calculated directly without explicitly calculating Gl and 
Gr,. In situations where also the SGF are needed, these 
can be obtained by using the relation 



Gl = -[Ko + ^l] , 
Gr = -[A-q + Er]-^ 



(46) 
(47) 



This can be derived by adding one layer to the left and 
one to the right surfaces respectivelyi^ In Appendix A 
we show that the SE calculated with eqs. (|44)) and ([iS]) 
indeed fulfill the above equation. Moreover with the use 
of eqs. and we can now regularize equation ([77]) 
also for the case where Tr is singular by writing it as 



9oo = ^^^0 - El - Er. 



(48) 



We have therefore a scheme where the SE are identified 
as the principal quantities, whereas the SGF and goo are 
derived from these. 

When we compare the method of reference [l| with 
the equations derived above, we notice that now it is not 
necessary to calculate the matrix goo and its inverse using 
eq. ((27)) in order to obtain the SE. This is not defined in 
the case of singular Ki and K^i, and therefore we expect 
the new method to be more stable and accurate. Also the 
problems caused close to band edges by the Van Hove 
singularities in goo are avoided. Moreover the method in 
reference [l[ relies on the calculation of the SGF in order 
to obtain the SE, whereas here the SGF is not needed. 
As we will show in section I VI I close to surface states the 
error in the SGF is much larger than the one for the 
SE, so that we also expect a large improvement in the 
accuracy for those particular states. 



IV. REDUCING THE CONDITION NUMBER 

OF A'l AND K-i 

The accuracy with which the SE are calculated de- 
pends on the accuracy involved in solving eq. ([6]), a 
quadratic eigenvalue problem extensively studied in the 
past Jiii^ However most solution methods have problems 
if Ki or K-i are close to being singular, or more gen- 
erally if their condition number k is large. In this case 
some of the complex eigenvalues tend to infinity and oth- 
ers to zero at the same time, and this results in a loss of 
accuracy in numerical computations. When calculating 
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FIG. 2: Schematic diagram of the basic algorithm described in section HITl and of the extended algorithm described in section 



Tr (7r) however the contributions from the states with 
Ini(fc„) — > oo (Im(fc„) —oo) are vanishingly small. 
It is therefore useful to limit the range of the eigenval- 
ues {e*'"'"} in such a way that the important eigenstates 
with small |Im(fc„)| and |Im(A:„)| can be calculated accu- 
rately, while losing precision for the less important eigen- 
states with large |Im(fc„)| and |Im(fc„)|. In this section 
we show how this can be achieved by decreasing k{Ki) 
and k{K_i). Here we assume that A'l K^-i, so that 
k{Ki) = k{K-i). Minor modifications are needed for the 
general case (see Appendix B). 

In order to obtain k{Ki) first a SVD of the matrix is 
performed 

Ki = USVl (49) 

U and V are unitary matrices, and 5" is a diagonal ma- 
trix, whose diagonal elements s„ are the singular val- 
ues. These are real and positive, and ordered so that 
Sn+i < Sn- If Smax IS the largest singular value, and Smin 
the smallest one, then the condition number is defined as 
k{Ki) = Smax/smin, with Ki singular if Smin is zero. 

We now replace S with an approximate S'svd, whose 
diagonal elements ssvD,n are 

1^ Ssvd Sn < Smax OSVD 

and accordingly Ki with A'l^svD — -S'svd The tol- 
erance parameter 6svb is a real positive number that 
determines the condition number of i^i,svD- 

We now present two possible choices for sgvo- The 
first is to set ssvd = 0, resulting in ATi^svD being sin- 
gular. We can then perform a unitary transformation 



in order to eliminate the degrees of freedom associated 
to ssvD.n = 0, and obtain an effective Ki matrix (ATJ^^) 
with reduced size for which k(A'°*) < (^svd- ^^^'^ second 
possibility is to set ssvd = Smax <^svd, so that by defini- 
tion we have k{Ki) — SgyY). The accuracy obtained with 
both strategies is similar, the advantage of using the first 
however is that the size of the matrices is reduced, so 
that for big systems the computation is much faster. In 
our implementation we use both methods together, first 
we reduce the size of the system by setting sgvD = 0, 
and then, if necessary, we further reduce the condition 
number for the effective system by limiting the smallest 
singular value. 



A. Reduction of system size 

Here we set all the M singular values s„ smaller than 
Smaxf^svD to zcro, SO that there are A^cft = N — M sin- 
gular values Sn with s„ > Smaxi^svD- The transforma- 
tions needed in order to obtain the right SE are now 
presented (the procedure for the left SE is analogous). 
We apply the unitary transformation K'^^, = U^Kzz'U, 
'?^R,« = U'^4'K,n, and we define K[ = WKi^syuU, K'_^^ = 
U^K^i^SYuU, Kq = WKqU. Since M singular values 
of A'l svD are zero the transformed matrices have the 
structure 
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where the dimensions of the new matrices are: A^eff x -^eff 
for Ki^c, K_i^c and A, N^e x M for Ki^^ and B,Mx N^s 
for /v_i^u and C, and M x M for D. Finally 0c, « is a 
column vector of dimension A'off , and 0u,„ is of dimension 
M. The transformed form of eq. (|4]) is 

(X^ + A-e»^" + A- le-'^-") 0j,,,„ = 0. (52) 

Due to the structure of K'_^ there are M solutions to this 
equation with e**"'" = and 4>c.n = 0. We therefore split 
up the right-going states into those with finite e*'^" ^ 
and those with e*'"'" = 0. For the first set, from eq. ([5^ . 
we obtain 

4>u,n = Fn4>c,n, (53) 

with 

F,, = -D-^ (JC-i,ue-"=" + C) . (54) 

The 0c, n are then solutions of an effective system with 
reduced size 

{Kf + Kfe"'- + Kf.e-^''") 0c,„ - 0, (55) 

where the effective matrices are 

K'^\ ^ /f_i,c - BD-^K^i,^, (56) 
Kf = A- BD-^C - Ki^nD-^K-i^n- 

We can now solve the quadratic eigenvalue problem (eq. 
([6])) for this effective system to get the set of A'^eff eigen- 
vectors Qc = ( 0c, 1 0c, 2 • ■ . 0c,Afetf ) and eigenvalues 
|gifc„| ^]^g right-going states. The M eigenvectors of 
the second set of solutions with e*'^'" = are given by 
0c, n = with a general 0u.n- The set of eigenvectors of 
the full AT' matrix therefore is 

with Qu = ( Ai0c,i A20c,2 •■• FAr^ff0c,Arrff ), and Qo is 
a general matrix of solution vectors for the states with 
gifeji _ From this we obtain the set of duals 

= ( -Qo%uQ-' Qo' ) • ^^^^ 

Using these results we can now calculate the transfer ma- 
trix of the transformed system 




where we have also used the fact that e'*^" = for the sec- 
ond set of solutions. We note that setting the M smallest 
singular values s„ to zero causes the last M columns of 
to be zero too. Moreover the explicit calculation of 
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FIG. 3: Absolute value of the diagonal elements of the 

transformed right SE for different values of (5svd,o. 

Qo is not needed in order to obtain T^. From this and 
eq. we obtain the right SE 

where 

= 5]e'^-"0c,„0j,„ (61) 

n=l 

is the SE of the effective system. 

The structure of shows that by applying this uni- 
tary transformation we have ordered the elements of the 
SE by absolute size, moving those columns (rows) with 
the smallest values to the right (bottom). By setting the 
smallest singular values of ATi to zero those columns and 
rows of the SE with small values have also been set to 
zero. This is illustrated in figure [31 where the absolute 
value of the diagonal elements of the transformed selfen- 
ergy |Sp | is shown for a (8,0) zigzag carbon nanotube at 
the Fermi energy Ap (sec section|V]for a detailed descrip- 
tion of the system). The are basically identical for 
different (5s vd up to i = TVoff, and indeed an increasing 
value of (5svD results in more diagonal elements of S'^ set 
to zero. We note that Nes is of similar size as N in figure 
[31 since the system is rather short along z and a small 
basis set is used (i.e. N is small). For large systems and 
rich basis sets the ratio N^s/N will decrease. The physi- 
cal interpretation of the zero columns and rows in the SE 
is that the M states with kn oo decay infinitely fast, 
so that the interaction of those states is limited to the 
site they are localized at. Finally the SE of the original 
system can be obtained by applying the inverse unitary 
transformation 

Sr = U^'^U^, (62) 

and in contrast to the matrix Sr is a dense N x N 
matrix. 

Note that in order to obtain the left SE we perform the 
unitary transformation K'^^, = V^Kzz'V, 0r „ = ^0r „, 
and then follow an analogous procedure. In this case 
however instead of the right-going states the left-going 
ones are used. 



8 



B. Limiting the smallest singular value 



V. ERROR ANALYSIS 



We can limit the lower bound of the singular values 
s„ by setting ssvd = SmaxfevD in eq. ((50)) . In this case 
the approximated K matrix is obtained by replacing Ki 
with Ki^svD- The error introduced is now of the order 
of Smax (^SVD- Ideally Smax (^SVD should be of the order 
of the machine numerical precision, so that the error is 
minimal. However sometimes increasing s,nax i^svd be- 
yond that value improves the results, therefore (5svd is 
left as a parameter to adjust depending on the material 
system investigated. This will be discussed extensively 
in the next section. 

A simpler but equally effective possibility for limiting 
the smallest singular value of a matrix is that of adding 
a small random perturbation . ^'^'^'^ Thus another strategy 
for reducing the condition number of Ki is that of replac- 
ing Ki with A'l^noisG = Ki+Wi^Wnoise), whcrC W{Wnoisc) 

is a matrix whose elements are random complex num- 
bers with an average absolute value \Wij\ ~ Wnoisc- In 
particular we choose the | Wij \ in such a way that both 
Re(Wij) and Im(Wij) are random numbers in the range 

[-WnoiscU'noisc]- Wc find that if tilnoisc = ^max (^SVD thc 

addition of noise usually gives results as accurate as those 
obtained with thc SVD procedure, but the calculation is 
faster since instead of performing a SVD we just perform 
a sum of the matrices. 



In figure [2] we present our final extended algorithm as 
it has been implemented in Smeagol. This now includes 
the following regularization procedure of Ki. First the 
size of Ki and hence of the whole problem is reduced 
by using the scheme described in section IIV Al with a 
tolerance parameter Jsvd = <5svd,i- This generates an 
effective matrix Kf^ whose condition number K{K'j^) is 
reduced by adding a small noise matrix VF(wnoisc)- Such 
a step is extremely fast and enhances considerably the 
numerical stability of thc calculation. In most cases the 
SE for the effective system can then be calculated and 
no further regularization steps are needed. However, in 
some cases the calculation of the SE still fails. This, for 
example, happens when the solution of eq. ^ for the 
effective system fails, or else when the calculated number 
of left-going states erroneously differs from the number of 
right-going states. In these critical situations wc further 
decrease k{K^^) by limiting the smallest singular value 
of Kf^ as described in section HV Bl with a tolerance pa- 
rameter ^svD = '5svD.2- The code automatically adjusts 
(^SVD,!, '5svD,2 and Wnoiso within a given range until the 
SE is calculated. In our test calculations for a number 
of different systems wc found no situation where such a 
scheme has failed. In contrast when thc standard algo- 
rithm of reference [l| is employed the number of failures 
was considerable. Note that our extended algorithm can 
also be used in conjunction with recursive methods for 
evaluating the SE.— i^ii^ii^ Also in this case it will decrease 
the computing time for large systems due to the reduced 
size of the effective K matrix. 



When recursive algorithms are used the accuracy of the 
SE is automatically known as it coincides with the con- 
vergence criterion. Poor convergence is found when the 
error can not be reduced below a given tolerance. Direct 
methods, as the one presented here, are in principle error 
free in the sense that when the solution is found, this is 
in principle exact. For this reason thc numerical errors 
arising from direct schemes usually are not estimated. In 
this section we perform this estimate and present a de- 
tailed error analysis for three different material systems. 

In order to estimate the numerical accuracy we use the 
recursive relations of eqs. (|46)) and (|47|) . written as 



Ki 



Vout 



(63) 



where E^^yj^j are calculated with our extended algo- 
rithm, and E'l^^^yp^j are obtained by evaluating the right- 
hand side term of the above equations. When the solu- 



tion is exact then E£"* 



E[" and E^"* 



we can define a measure of the error as 



EJ^\ Therefore 



^{L/R} 



^{L/R} 



(64) 



where 



stands for the max norm,— thc corre- 



sponding relative error is A^^ 



Ae/||e 



{L/R} I 



The accuracy criterion used in the extended algorithm 
is the following. We first set (5svd,i, Wnoiso and eventu- 
ally (5svD.2 and compute As^i-. This should be lower than 
a target accuracy A^'^^. If this is not the case then the 
SE will be recalculated with a different set of tolerance 
parameters, until A^.r reaches the desired accuracy. If 
this condition is never achieved the final SE is thc one 
with to the smallest As^r . 

We now calculate the SE for different variations of the 
method, chosen in order to highlight the problems aris- 
ing from Ki and K-i and to show the difference between 
the basic method of reference [l| and thc extensions pre- 
sented here. There arc two main differences between the 
two methods. The first is that here we solve eq. ([6]) with- 
out inverting Ki, whereas in reference [l[ K^^ is used 
to solve the inverse band-structure relation k = k{E). 
Clearly this second choice is less accurate if Ki is close 
to singular. However it is much faster computationally, 
so that it might be of advantage for big systems. The 
second difference is that here it is not necessary to calcu- 
late goo via eq. (j27p . so that one does not need to invert 
Tr, and Tr. 

In order to investigate the effect of these two aspects 
independently, we have calculated the SE using the fol- 
lowing four methods. In method 1 we use the algorithm 
presented in this work. In particular we use eq. ([6]) 
to solve the quadratic eigenvalue problem and eqs. (Pil - 
HS]) to obtain thc SE (for thc right SE wc actually use 
a different form of eq. sec Appendix C). Method 2 





FIG. 4: Unit cells of the three systems investigated in this 
work: (a) (8,0) zigzag carbon nanotube, (b) bcc Fe oriented 
along the (100) direction, (c) fee Au oriented along the (111) 
direction. The black arrow indicates the direction of the 
stacking 2:, i.e. the direction of the transport. 



is essentially the same, with the only difference that in- 
stead of solving eq. (O we use the eigenvalue method of 
reference In method 3 we solve eq. ([5]), but we use 
eq. ((42l) to calculate the SGF, with goo obtained from eq. 
([j?]) . Finally method 4 is the algorithm of reference 

In order to obtain a statistically significant average of 
the errors, we plot a histogram of the calculated errors for 
both El and Sr for a large energy range. Here we use 
the absolute error, since it can readily be compared to 
the energy scale of the problem. Note that although the 
relative error might be small, the absolute error can be 
very large if | |l]{L/pj j 1 ^ 1 Ry. Furthermore in order 
to keep the analysis simple in all the calculations of this 
section we do not reduce the system size nor do we add 
noise (wnoisc = 0). We regularize Ki and K-i by using 
ssvD = Sinax ^SVD in eq. ([50)1 . Since the error depends on 
the chosen Ssvd, here we calculate As for a set of <5svd 
in the range [0, 10-^3, IO-22, . . . , lO-^, IQ-^]. We then 
present the smallest As found for i5svd taken in that 
range. This is the smallest possible error achievable with 
a given method and allows us to extract informations on 
the range of optimal SVD values for a given method. 

As first example a (8,0) zigzag carbon nanotube^ is 
presented (the unit cell is shown in figure H^a)). The 
length of the periodic unit cell is 4.26 A along the nan- 
otube, with 32 carbon atoms in the unit cell. The LDA 
approximation (no spin-polarization) is used for the ex- 
change correlation potential. We consider 2s and 2p or- 
bitals for carbon with double-^ and a cutoff radius Tc for 
the first C of = 5 Bohr. Higher C are constructed with 
the split-norm scheme with a split-norm of 15%fi^ The 
real space mesh cutoff is 200 Ry. The matrices Hq, Hi, 
Sq and Si are extracted from a ground state DFT calcu- 
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FIG. 5: Histogram of the errors in the calculation of the self- 
energy Ae for three different systems, (a) (8,0) zigzag carbon 
nanotube, (b) bcc Fe, (c) fee Au. A'^ is the number of times a 
given error Ae occurs (not normalized). 
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FIG. 6: Histogram of 5s vd giving the smallest error in the 
self-energy, (a) (8,0) zigzag carbon nanotube, (b) bcc Fe, (c) 
fee Au. A'^ is the number of times a given Ssvo generates the 
smallest error (not normalized). 



lation for an infinite periodic nanotube. We calculate the 
SE for the semi-infinite nanotube at 1024 energy points 
in a range of ±5 eV around the Fermi energy. 

Figure [DJa) shows the histogram of the errors in the 
SE, where is the number of times a given error As 
appears. In general the figure shows that for this system 
the average error increases when going from method 1 to 
method 2 and method 3, and finally to method 4. The 
error obtained with method 1 is on average about 6 orders 
of magnitude smaller than the one obtained with method 
4. The main reason behind this dramatically improved 
accuracy is that method 1 does not involve any steps 
where a singular Ki leads to divergencies. Method 4 on 
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the other hand is strongly dependent on the condition 
number of Ki , since it necessitates to invert Ki and Tr 
(or Tr). Methods 2 and 3 are on average about one order 
of magnitude more precise than method 4. Since they 
both still involve one of the two inversions the difference 
is however not large. 

Figure [6ja) shows the histogram of the optimum 5svd 
used for the calculations of the SE. Here we plot the num- 
ber of times N a particular Ssvu has given the smallest 
error in the set of calculations. A larger optimal value 
for (5svD indicates a stronger dependence of the computa- 
tional scheme on ^(i^i). For method 1 the range of used 
<5svD is smaller than 10"-'^^. If we force Ssvb to be zero we 
get almost the same level of accuracy as shown in figure 
[SJa), which confirms that the accuracy of for method 1 
depends little on k{Ki) for this system. However also for 
this method there is a set of energies (a few percent of the 
total number) where the solution of eq. Q fails if Ssvd 
is too small. The optimal Ssvd for the other methods is 
orders of magnitude larger than that of method 1 , and it 
is never smaller than 10~^. The absolute error induced 
by replacing Ki by Ki^syb is of the order of (Ssvd Smax- 
Usually Smax is of the order of 1 Ry, so that the error is 
of the order of 6syu Ry- Therefore since in methods 2 
to 4 a large value of i5svd is needed in order to improve 
^(ii'i^svD), also the resulting error is large. 

The second example is bee Fe (figure Hl[b)), oriented 
along the (100) direction. The lattice parameters are 
the same as in reference [1^. There are 4 Fe atoms in 
the unit cell. We apply periodic boundary conditions in 
the direction perpendicular to the stacking, so that these 
correspond to 4 Fe planes. The length of the cell along 
the stacking direction is 5.732 A. A doublc-C s (rc=5.6 
Bohr), single-C p (rc=5.6 Bohr) and single-C d (rc=5.2 
Bohr) basis is used. The real space mesh cutoff is 600 Ry, 
and the DFT calculation is converged for 7x7 fc-points in 
the Brillouin zone orthogonal to the stacking. The SE 
have been calculated for the converged DFT calculation 
at 32 different energies in a range of ±1 cV around the 
Fermi energy, and for 10,000 /c-points in the 2D Brillouin 
zone perpendicular to the stacking direction. For each 
fc-point there is a different set of matrices Kq, Ki and 
K-i, so that for each fc-point there is a different SE. The 
histogram for the error of the calculated self-energy As 
is shown in figure[n{b), and the histogram for the optimal 
SsvB in figure [SJ^b). The general behavior is similar to 
the one found for the carbon nanotube. We note that, al- 
though for the vast majority of the calculations the error 
in the SE is small, there is a long tail in the histograms 
of figure [5jb) indicating the presence of a small number 
of large errors. This is present for all the methods, with 
a maximum error of 10~^ Ry for method 1, and 100 Ry 
for method 4. Closer inspection shows that the reason for 
the increase of the error for certain energies and fc-points 
is caused by a divergence in | |S{l,r} | l^j^^' This will be 
illustrated in more detail in the next section. 

Finally we consider fee Au (figure [3)Jc)), with the stack- 
ing along the (111) direction. The miit cell consists of 



three planes of nine gold atoms each. These are the typ- 
ical leads used for the calculations of the transmission 
properties of molecules attached to goldi ^^i^'^'^^i^^ We use 
double-C s (rc=6.0 Bohr) and single-^ d (rc=5.5 Bohr) 
and four fc-points in the Brillouin zone perpendicular to 
the stacking. The mesh cutoff is 400 Ry. The SE have 
been calculated for 418 energy points, from about 15 eV 
below to about 10 cV above the Fermi energy. The gen- 
eral behavior (figures [Sljc) and[6l[c)) is again similar to 
that of the previous examples. Also here the error for 
method 1 is about 6 orders of magnitude smaller than 
that of method 4, with method 2 and 3 giving some 
marginal improvement. 

Our results show that the new scheme in general al- 
lows the calculation of the SE with high accuracy. The 
main advantage of method 1 is rooted in the possibility 
of using a much smaller (Ssvd • For big systems sometimes 
one might prefer to use method 2, since it is considerably 
faster than method 1 and gives the second best accuracy. 
In this case we first calculate the SE with method 2 and 
check the error. Only for those energy points where the 
error is above some maximum value (of the order of 10~^ 
Ry for example) the calculation is repeated with method 
1 to improve the accuracy. Finally the results show that 
for all methods the SVD transformation of A'l is nec- 
essary, although for method 1 it is needed only a few 
percent of the times. For big systems, in particular if the 
unit cell is elongated along the stacking direction, or if 
a rich basis set is used, k{Ki) will generally increase as 
there will be some singular values of Ki going to zero. 
In these cases also method 1 will require a SVD trans- 
formation for most energies. The range of (Ssvd should 
however be similar to the one shown in figure [6l so that 
also the error in the SE should be of the same order of 
magnitude. We also note that in order to keep the analy- 
sis simpler here we have not used the reduction of system 
size described in section HV A[ for such large systems it 
is however crucial in order to decrease the computational 
effort and regularize Ki at the same time. 



VI. SURFACE STATES 

The center of the error distribution for method 1 (fig- 
urc[5]) is located at small As, usually smaller than 10~^^ 
Ry. However the histogram has also a tail reaching up to 
very large errors. These are found only at some critical 
energies as demonstrated in figure [TlJ a), where we show 
As for the carbon nanotube calculated over 1024 energy 
points in a range of 2 cV around the Fermi energy. The 
average error is of the order of 10~^^ Ry, but at ener- 
gies around -0.8 eV and -0.34 eV the error drastically 
increases. Indeed a finer energy mesh at these points 
suggests a divergence. The origin of the large errors at 
particular energies can be investigated by looking at the 
eigenvalues gh,i of the SGF Gl- In figure [Tljb) the largest 
and the smallest absolute value for the eigenvalues, re- 
spectively 5L,max and 5L,min, are plotted as function of 
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energy (^L.min < \gL,i\ < 5L,max)- It can be seen that 
SL,max diverges close to the energies where the error in- 
creases, i.e. we can associate large errors in Gl with a 
divergence in its spectrum. Since El is calculated from 
eq. ([44]) the only possible origin for the divergence is in 
the norm of some of the (j)-R,^n- As these are obtained by 
inverting the matrix Q = ( 0r,i 0r^2 • ■ • 4>R,n ) (cq. 
(fTO)) ). one deduces that the set of vectors {0r.„} is not 
linearly independent. For these energies k{Q) oo. We 
therefore can simply check the magnitude of k{Q) to de- 
termine whether there is a divergence of the SE close to 
a particular energy. 

Physically the divergence of the SE translates into the 
presence of a surface state at that particular energyi^ 
Consider the spectral representation of Gl 
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E + i5-E, 
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(65) 



where £"„ are the eigenvalues and ipn arc the right 
eigenvectors of the effective surface Hamiltonian matrix 
Hq — Sl with overlap 5*0, and tpn are the left eigen- 
vectors of the same Hamiltonian. A localized surface 
state is found when there is a real eigenvalue En{E) 
at En{E) = E (or more generally if Im(i?„(£')) is very 
small) . 

From the recursive relation (I46p one can deduce that 
for an infinite eigenvalue there is also a corresponding 
vanishing eigenvalue. Therefore in figure [7{b) for ener- 
gies where (/L,max — ^ ^ wc have also </L.min — ^ 

0. Close 

to the singularity we can therefore expand the two eigen- 
values as gL,max OC and 5L,min oc E + iS - En- 

For E = En the largest eigenvalue in eq. ([65]) is then 
equal to 6~^, and the smallest is equal to 5. To avoid 
divergence therefore the magnitude of the Gl eigenval- 
ues can be bounded to a finite value by introducing 
a small imaginary part to the energy for energies in the 
vicinity of a surface state. 

Another possibility for limiting the size of (7L,max is to 
bound the singular values of Q from below in the same 
way as it is done for Ki (section llV Bp . This essentially 
imposes the 4>B..n to be linearly independent from each 
other. However, with this scheme it is not possible to 
conserve the Green's function causality, so that the SGF 
might have eigenvalues lying on the positive imaginary 
axis. Moreover we loose control over the accuracy of the 
computed SGF and SE. Both these problems are avoided 
when using a finite S. 

We now investigate the DOS and transport properties 
of a system when the finite imaginary part S (broadening) 
is added to the energy. We consider as an example the 
carbon nanotube of figure 21 In figure [5Ja) the onsite 
surface DOS po as defined in eq. (|551) is shown for 6 = 
Ry, 6 = 10-6 Ry, S = 10^5 Ry and b = 10"^ Ry 
For (5 = the surface DOS vanishes for energies between 
-0.37 cV to +0.45 cV, indicating the presence of a gap 
around the Fermi energy. Note that there are no Van 
Hove singularities in po, since we never divide by the 
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FIG. 7: Error analysis for the carbon nanotube of figure^ (a) 
absolute error Ae of the self-energy as function of the energy 
-E, (b) maximum ((jL.max) and minimum (gL.min) eigenvalues 
ofGL. 
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FIG. 8: Density of states and transmission coefficient for the 
carbon nanotube of figure [l] (a) Density of states of the 
surface layer po as function of energy _E, calculated for dif- 
ferent broadenings 5. The inset is a zoom at energies around 
-0.34 eV. (b) Transmission coefficient T for different values of 
5. 



group velocity when calculating the SGF. For finite (5 and 
energies away from the band gap, the DOS is essentially 
identical to that calculated for (5 = 0, however inside 
the gap pq does not vanish but saturates to a small value 
proportional to . Moreover whereas the surface states 
are not visible for (5 = 0, they appear in the DOS for finite 
(5, with a full width at half maximum (FWHM) equal to 
2(5. 

We then move to the transport by calculating the 
transmission coefficient'^ T(E) for a carbon nanotube at- 
tached to semi-infinite leads made from an identical car- 
bon nanotube. Since this is a periodic system T{E) must 
equal the number of open channels, so that it can only 
have integer values. This is indeed the case for 6 = 
(figure [H{b)). For finite 6s the transmission coefficient is 
only approximately an integer, especially inside the en- 
ergy gap region where the finite surface DOS introduced 
by 6 leads to a non zero transmission. The transmission 
in the gap is proportional to 6^ (note that the scale is 
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logarithmic) , since on both sides of the scattering region 
the artificial surface DOS is proportional to i5. In this 
region of small transmission therefore the results might 
change by orders of magnitude depending on the value 
of 6. For all values of 5 however wc find no contribu- 
tions to the transmission coming from the surface state, 
indicating that these do not carry current. These results 
show that adding a finite value S to the energy has little 
effect on the actual transmission if this is large. However 
when the transmission is small, as in the case of tunnel 
junctions, the finite S introduces an additional contribu- 
tion to the conduction that might arbitrarily affect the 
results. It is thus imperative for those systems to identify 
surface states and use the imaginary S only in a narrow 
energy interval around them. 

Finally we can give an estimate of the relative accuracy 
^s,i{S) ~ As/llSll^ax ^^'^ energy corresponding to 
the surface state. As discussed before the origin of the 
error is the inversion of Q needed to calculate the duals. 
The relative error introduced by the inversion of Q is 
proportional to K((5)fi^iiii^i^i^ Close to a surface state 
the smallest singular value is of the order of 6, so that 
k(Q) cx S~^. As this is the dominant source of error in 
the calculation of the SE close to a surface state, we can 
approximate the relative error as 



ci S 



(66) 



where ci is a constant that depends on the machine pre- 
cision and on the details of the algorithm. The label "in" 
explicitly indicates that this is the error in the SE calcu- 
lated with the extended algorithm in eq- (|63p ). 
The absolute error A'^P is equal to the relative error times 
ll^llniiix' which is itself proportional to S~^, so that we 
get A^^ cx S-^. 

When using eq. to estimate the error in the SE 
we introduce an additional error due to the inversion in- 
volved in obtaining Gl • The largest singular value of Gl 
is proportional to S^^, and the smallest one is propor- 
tional to S, so that the relative error introduced by the 
inversion is proportional to k{G£^) = k(Gl) oc 5"^. For 
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where C2 is again a constant. Since the errors are random 
the total estimated error can be approximated by adding 
the contributions from the two inversions 
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(68) 



As,r is therefore a good estimate for the true error A^'j. 
if A|;"J is small. Close to surface states however A™J » 
Ag ^, so that As^r largely overestimates the true error. 

To verify these estimates numerically we present a 
scheme for calculating A^ j. and A|,"* independently. For 
each SE we perform a second calculation where we add 
a small amount of noise to the input matrices Ko,Ki, 
and -ftT-i, so that we obtain the self-energy Sl, noise for 
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FIG. 9: (a) Relative error of the self-energy (A^,r represents 
the true error), (b) condition numbers of Q and Gl, as a 
function of the broadening 5 for the carbon nanotube of figure 
|4] calculated at the surface state energy. 



a slightly perturbed system. The noise is added as a 
random relative perturbation of each element of the ma- 
trices. As we decrease the magnitude of the noise the 
difference between El and Sl, noise is reduced until it be- 
comes constant for noise smaller than a critical value. 
In this range of minimum noise even if the difference 
in the input matrices decreases, the difference in the 
output matrices is constant, it therefore corresponds to 
the error in the calculation. As one might expect we 
find that this critical value of noise is of the same or- 
der of magnitude as the numerical accuracy used (ap- 
proximately 10"^^ in our calculations). We can there- 



fore obtain A!; 



7 



E;" and 

^ I I max 



A out 



— II^L 



_ yout I 
L, noise | 



"'L, noise 

/I IE?"* 1 1 , with the mag- 

'II L II max ' o 

nitude of the noise equal to the critical value. 

We have calculated the maximum error for a set of 128 
energy points located within 10^^^ Ry around the energy 
of the surface state at -0.34 eV for different values of 5. 
The result is shown in figure [^U a). Indeed for small S 



A IE 

eq. 



follows eq. 
67]) with C2 



Ml) with ci « 10-^^ Ry, A' 
i 10-19 Ry2, and (As,r)^ « (A 



follows 

(A|;"J:)^. In figurcinUb) the condition numbers k{Q) and 
k(Gl) are shown, confirming k{Q) oc S^^ and k(Gl) oc 
S^^. This demonstrates that close to surface states As r 
is mainly caused by the calculation of Gl. Thus As r 
largely overestimates the real error A^^, which even for 
S = IQ-i" Ry has an acceptable size of A^\. ~ 10~^. 

Since ci and C2 are generally system dependent, in 
practical calculations we use a value of 6 ranging between 
10"^ Ry and 10"® Ry for energies in the vicinity of sur- 
face states, mainly in order to limit the absolute error. 
Moreover S is added in an energy range corresponding 
approximately to the FWHM of the imaginary part of 
{E — En + iS)~^, which is equal to 26. Although this 
range is only of the order of 10"'' — 10"® Ry, in practical 
calculations where both energy and fc-point sampling is 
fine the number of times when this prescription is applied 
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can be rather large (see figure [5]). 

The above analysis confirms that close to surface states 
also direct methods have the same accuracy problems of 
recursive methods. This fact is usually ignored in the 
literature,— i^iiSi^ where it is assumed that the accuracy 
is constant for a given algorithm. Here we show that the 
accuracy of a method is solely determined by the value 
of Ci , which, as indicated in section |Vl can vary over 
many orders of magnitude. Our analysis also shows that 
methods requiring the explicit calculation of Gl from its 
inverse are much less accurate close to surface states than 
those calculating El directly. 



VII. CONCLUSIONS 

By extending the scheme proposed in reference [l| we 
have presented a different but equivalent form for calcu- 
lating the Green's functions of an infinite quasi-lD sys- 
tem, as well as the SGF and SE for the semi-infinite sys- 
tem. We have then constructed an extended algorithm 
containing also the necessary steps to regularize the ill 
conditioned hopping matrices. This is found to be cru- 
cial in order to obtain a numerically stable algorithm. 
By applying a unitary transformation based on a SVD 
we remove the rapidly decaying states and calculate the 
SE for an effective system with reduced size. We further 
decrease the condition number of the hopping matrices 
by adding a small random perturbation and by limiting 
the smallest singular value. 

We have performed a detailed error analysis on the 
numerical calculation of the SE, showing that if the algo- 
rithm does not involve an inversion of the hopping ma- 
trices Ki (or K-i) high accuracy is obtained. We also 
find that the error is not constant as function of energy. 
It is shown that an increase of accuracy is needed espe- 
cially close to energies where the SE and SGF diverge, 
which corresponds to the presence of surface states in the 
semi-infinite system. At these energies we improved the 
accuracy by adding a small imaginary part to the energy. 
We have shown that this procedure affects the transport 
properties little in the high transmission limit. How- 
ever, for low transmission this adds some spurious sur- 
face density of states contributing significantly to the to- 
tal transmission. The transport can therefore be strongly 
affected, so that the imaginary part should be added only 
in a small energy range around the poles and it should 
be as small as possible. 

Our final algorithm is therefore highly numerically sta- 
ble and extremely accurate. Most importantly errors and 
accuracy can be closely monitored. We believe that this 
is an ideal algorithm to be used with ab initio transport 
schemes, where the condition of the Hamiltonian and its 
sparsity is controlled by the convergence of the electronic 
structure and therefore cannot be fixed a priori. 
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APPENDIX A: VERIFICATION OF THE 
RECURSIVE RELATION FOR THE SGF 

Here wc demonstrate that El calculated using eq. (144]) 
indeed fulfills the recursive relation for Gl of eq. (|46p . 
Insert eqs. ([^S]) and (jH]) into eq. and take the inverse 
to obtain 

Ko + K^iTn + KiT^' = 0. (69) 
Using the definition of the matrix Tr (eq. (|24|) ) we write 

N 



^ ^R,n4,„ = 0. (70) 



This equation corresponds to the defining equation for 
the (^R,n and is therefore fulfilled by definition. The same 
is therefore true for eq. Eq. (|T7|) for Gr. can be 

demonstrated similarly. 



APPENDIX B: REGULARIZATION OF Ki AND 

K-i FOR Kl ^ K-i 

In section HV Al we assume that A'l = K^i in order to 
write the transformed matrices K[ and K'_i in form of 

eq. ([5T|). If kI ^ K_i the same can be done by perform- 
ing a generalized SVD of the Hamiltonian and overlap 
matrices as described in reference Here we present a 
different approach, based on two standard SVD transfor- 
mations, one for Ki and one for K^_i 



(71) 



Here J7i,C/_i,Vi and V-i are unitary matrices, S'a and 
Sh are diagonal matrices with the singular values on the 
diagonal. In general there are Mi singular values of Ki 
smaller than (5svDSa,max, and M_i singular values of K-i 
smaller than i5svDSb,max, with Sa.max and Sb.max being 
respectively the largest singular value of Ki and 
If M = min(Afi, M_i), wc obtain Ki^syu by setting the 
smallest M singular values of Ki to zero. In the same way 
we obtain K^i^gVB by setting the smallest M singular 
values of A'_i to zero. A transformation 



K[ = C//ifi,svDC/_i, 
K'_i = C///f-i,svDC/-i 



(72) 
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brings both K[ and K'_^ to the form of eq. ((5T|) . All the 
results of section IIV Al are then valid also for k\ ^ K- 1 . 

If the Hamiltonian and overlap matrices are real and 
Hermitian, but the energy is complex, then Ki = A'l* . 
By using eq. (ffTj) . and the fact that Sa and Sh arc real, 
we obtain = Sh, so that M = Mi = If the 

Hamiltonian and overlap matrices are Hermitian but not 
real, then in general Sh- However in all the calcu- 

lations performed the difference between Sa and Sh was 
very small, so that in practice we always had Mi = M2. 

In section HV BI we limit the singular values of Ki from 
below without reducing the size of the system. If A'J ^ 
K_i we simply apply the transformations described in 
section ITVB I to both Ki and K-i independently. 



APPENDIX C: QUADRATIC EIGENVALUE 
PROBLEM FOR THE RIGHT-GOING STATES 

We find that in the solution of eq. ([6|) the numerical ac- 
curacy for those eigenvalues with \e^''" \ > 1 (Im(/c„) < 0) 



is better than for those with je*'^"] < 1 (Im(fc„) > 0), 
especially when \kn\ ^ 1- For El we only need the left- 
going states, for which eq. ^ gives the better accuracy. 
For Er the right-going states are needed. In this case, 
in order to increase the accuracy for the right decaying 
states (Im(A:„) > 0), instead of eq. ([6]) we solve the equiv- 
alent equation 

with 

'J>H. = (74) 

The eigenvalues of the states with Im(fc„) > now have 
an absolute value larger than one and therefore a higher 
accuracy. 
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